*****************************************************************
*** Replication file for "Long-Term Effects in Models with Temporal Dependence", Political Analysis, Laron K. Williams
***
*** NOTE: Any questions, please email williamslaro@missouri.edu
*****************************************************************

*** Set up the working directory containing all the files
*cd "FOLDER LOCATION"

set seed 648

*****************************************************************
*** Monte Carlo experiments
*****************************************************************

* Produce the data with "Monte Carlo Experiments.R" (calling subprograms in "Monte Carlo Programs.R")
* Analyze the data (producing Table 1 and others in Appendix) with "Analyze Monte Carlo Experiments.dta"

*****************************************************************
*** Long-term effects: Negative duration dependence
*****************************************************************
use "Negative.dta", clear

btscs y tim group, g(t) nspline(3)
replace t = (t + 1)
gen t2 = t^2
gen t3 = t^3

set seed 648
estsimp logit y x t t2 t3

************************ All Types of Effects *****************************
* Rows: x (1-time), x (temporary), x (permanent), t0_0, t0_1, t1_0, t1_1
* Columns: Baseline, counterfactual, then simulations
mat S = (0 , .5,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0 \ /*
*/  	 0 , .5, .5, .5, .5,  0,  0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0 \ /*
*/		 0 , .5, .5, .5, .5, .5, .5, .5, .5, .5, .5, .5,.5,.5,.5,.5,.5 \ /*
*/		 13, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28 \ /*
*/		 13, 13,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14 \ /*
*/		 13, 13,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14 \ /*
*/		 13, 13,  0,  1,  2,  3,  4,  0,  1,  2,  3,  4,  5, 6, 7, 8, 9)

mat SIMS = S[.,3...]
local simno = colsof(SIMS)


tempname lteneg
postfile `lteneg' xaxis time0_0 time0_1 time1_0 time1_1 /*
*/	x1 pr_s1_0_0 pr_s1_0_0_lo pr_s1_0_0_hi pr_s1_0_1 pr_s1_0_1_lo pr_s1_0_1_hi pr_s1_1_0 pr_s1_1_0_lo pr_s1_1_0_hi pr_s1_1_1 pr_s1_1_1_lo pr_s1_1_1_hi lte_s1 lte_s1_lo lte_s1_hi /*
*/	x2 pr_s2_0_0 pr_s2_0_0_lo pr_s2_0_0_hi pr_s2_0_1 pr_s2_0_1_lo pr_s2_0_1_hi pr_s2_1_0 pr_s2_1_0_lo pr_s2_1_0_hi pr_s2_1_1 pr_s2_1_1_lo pr_s2_1_1_hi lte_s2 lte_s2_lo lte_s2_hi /*
*/	x3 pr_s3_0_0 pr_s3_0_0_lo pr_s3_0_0_hi pr_s3_0_1 pr_s3_0_1_lo pr_s3_0_1_hi pr_s3_1_0 pr_s3_1_0_lo pr_s3_1_0_hi pr_s3_1_1 pr_s3_1_1_lo pr_s3_1_1_hi lte_s3 lte_s3_lo lte_s3_hi /*
*/	using "LTE Example--Negative.dta", replace


qui foreach c of numlist 1(1)`simno' {
	local x1 = S[1,`c']
	local x2 = S[2,`c']
	local x3 = S[3,`c']
	local t0_0 = S[4,`c']
	local t0_1 = S[5,`c']
	local t1_0 = S[6,`c']
	local t1_1 = S[7,`c']
	
	*** Y_{t-1} == 0 & Y_{t-1} == 1	
	foreach s of numlist 1(1)3 {
		foreach b of numlist 0 1 {
			foreach o of numlist 0 1 {
				tempvar prob_s`s'_`o'_`b'
				setx x `x`s'' t `t`o'_`b'' t2 `t`o'_`b''^2 t3 `t`o'_`b''^3
				simqi, prval(1) genpr(`prob_s`s'_`o'_`b'')
				sum `prob_s`s'_`o'_`b'', meanonly
				local pr_s`s'_`o'_`b' = r(mean)
				_pctile `prob_s`s'_`o'_`b'', perc(2.5 97.5)
				local pr_s`s'_`o'_`b'_lo = r(r1)
				local pr_s`s'_`o'_`b'_hi = r(r2)
			}
		}
	
		tempvar lteff_s`s'
		gen `lteff_s`s'' = `prob_s`s'_0_1' - `prob_s`s'_0_0'
		sum `lteff_s`s'', meanonly
		local lte_s`s' = r(mean)
		_pctile `lteff_s`s'', perc(2.5 97.5)
		local lte_s`s'_lo = r(r1)
		local lte_s`s'_hi = r(r2)
	}
	
	local xaxis = `c' - 2
	
	post `lteneg' (`xaxis') (`t0_0') (`t0_1') (`t1_0') (`t1_1') /*
*/	(`x1') (`pr_s1_0_0') (`pr_s1_0_0_lo') (`pr_s1_0_0_hi') (`pr_s1_0_1') (`pr_s1_0_1_lo') (`pr_s1_0_1_hi') (`pr_s1_1_0') (`pr_s1_1_0_lo') (`pr_s1_1_0_hi') (`pr_s1_1_1') (`pr_s1_1_1_lo') (`pr_s1_1_1_hi') (`lte_s1') (`lte_s1_lo') (`lte_s1_hi') /*
*/	(`x2') (`pr_s2_0_0') (`pr_s2_0_0_lo') (`pr_s2_0_0_hi') (`pr_s2_0_1') (`pr_s2_0_1_lo') (`pr_s2_0_1_hi') (`pr_s2_1_0') (`pr_s2_1_0_lo') (`pr_s2_1_0_hi') (`pr_s2_1_1') (`pr_s2_1_1_lo') (`pr_s2_1_1_hi') (`lte_s2') (`lte_s2_lo') (`lte_s2_hi') /*
*/	(`x3') (`pr_s3_0_0') (`pr_s3_0_0_lo') (`pr_s3_0_0_hi') (`pr_s3_0_1') (`pr_s3_0_1_lo') (`pr_s3_0_1_hi') (`pr_s3_1_0') (`pr_s3_1_0_lo') (`pr_s3_1_0_hi') (`pr_s3_1_1') (`pr_s3_1_1_lo') (`pr_s3_1_1_hi') (`lte_s3') (`lte_s3_lo') (`lte_s3_hi')
}

postclose `lteneg'



*****************************************************************
*** Long-term effects: Negative duration dependence w/ Positive NPH
*****************************************************************
*** Modeled without NPH
use "Negative with Positive NPH.dta", clear

btscs y tim group, g(t) nspline(3)
replace t = (t + 1)
gen t2 = t^2
gen t3 = t^3
gen x_t = x*t

set seed 648
estsimp logit y x t

************************ Permanent Long-Term Effects *****************************
* Rows: x (permanent), t0_0, t0_1
* Columns: Baseline, counterfactual, then simulations
mat S = (0 , .5, .5, .5, .5, .5, .5, .5, .5, .5, .5, .5,.5,.5,.5,.5,.5,.5,.5,.5,.5,.5 \ /*
*/		 13, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33 \ /*
*/		 13, 13,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19)

mat SIMS = S[.,3...]
local simno = colsof(SIMS)

tempname lteneg_no
tempfile no_nph
postfile `lteneg_no' sim xaxis time0_0 time0_1 /*
*/	x1 pr_0_0_n pr_0_1_n lte_n /*
*/	using `no_nph', replace

qui foreach c of numlist 1(1)`simno' {
	local x1 = S[1,`c']
	local t0_0 = S[2,`c']
	local t0_1 = S[3,`c']

	local xaxis = `c' - 2
	
	*** Y_{t-1} == 0 & Y_{t-1} == 1	
	foreach o of numlist 0 1 {
		tempvar prob_`o'
		setx x `x1' t `t0_`o''
		simqi, prval(1) genpr(`prob_`o'')
	}
	local b = 1
	while `b' <= 1000 {
		local p0 = `prob_0'[`b']
		local p1 = `prob_1'[`b']
		local lte = `p1' - `p0'	
		post `lteneg_no' (`b') (`xaxis') (`t0_0') (`t0_1') (`x1') (`p0') (`p1')	(`lte')	
		local b = `b' + 1 
	}
}

postclose `lteneg_no'

*** Modeled with NPH
use "Negative with Positive NPH.dta", clear

btscs y tim group, g(t) nspline(3)
replace t = (t + 1)
gen t2 = t^2
gen t3 = t^3
gen x_t = x*t

set seed 648
estsimp logit y x t x_t

************************ Permanent Long-Term Effects *****************************
* Rows: x (permanent), t0_0, t0_1
* Columns: Baseline, counterfactual, then simulations
mat S = (0 , .5, .5, .5, .5, .5, .5, .5, .5, .5, .5, .5,.5,.5,.5,.5,.5,.5,.5,.5,.5,.5 \ /*
*/		 13, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33 \ /*
*/		 13, 13,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19)

mat SIMS = S[.,3...]
local simno = colsof(SIMS)

tempname lteneg_yes
tempname yes_nph
postfile `lteneg_yes' sim xaxis time0_0 time0_1 /*
*/	x1 pr_0_0_y pr_0_1_y lte_y /*
*/	using `yes_nph', replace

qui foreach c of numlist 1(1)`simno' {
	local x1 = S[1,`c']
	local t0_0 = S[2,`c']
	local t0_1 = S[3,`c']

	local xaxis = `c' - 2
	
	*** Y_{t-1} == 0 & Y_{t-1} == 1	
	foreach o of numlist 0 1 {
		tempvar prob_`o'
		setx x `x1' t `t0_`o'' x_t `x1'*`t0_`o''
		simqi, prval(1) genpr(`prob_`o'')
	}
	local b = 1
	while `b' <= 1000 {
		local p0 = `prob_0'[`b']
		local p1 = `prob_1'[`b']
		local lte = `p1' - `p0'
		post `lteneg_yes' (`b') (`xaxis') (`t0_0') (`t0_1') (`x1') (`p0') (`p1') (`lte')		
		local b = `b' + 1 
	}
}

postclose `lteneg_yes'
	
use `no_nph', clear

sort xaxis sim
merge xaxis sim using `yes_nph', sort

gen diff_lte = lte_y - lte_n
gen diff_0_0 = pr_0_0_y - pr_0_0_n
gen diff_0_1 = pr_0_1_y - pr_0_1_n
	
collapse (mean) time0_0 time0_1 x1 pr_0_0_y_mn = pr_0_0_y pr_0_1_y_mn = pr_0_1_y pr_0_0_n_mn = pr_0_0_n pr_0_1_n_mn = pr_0_1_n diff_0_0_mn = diff_0_0 diff_0_1_mn = diff_0_1 diff_lte_mn = diff_lte /*
*/	(p3) pr_0_0_y_lo = pr_0_0_y pr_0_1_y_lo = pr_0_1_y pr_0_0_n_lo = pr_0_0_n pr_0_1_n_lo = pr_0_1_n diff_0_0_lo = diff_0_0 diff_0_1_lo = diff_0_1 diff_lte_lo = diff_lte /*
*/	(p97) pr_0_0_y_hi = pr_0_0_y pr_0_1_y_hi = pr_0_1_y pr_0_0_n_hi = pr_0_0_n pr_0_1_n_hi = pr_0_1_n diff_0_0_hi = diff_0_0 diff_0_1_hi = diff_0_1 diff_lte_hi = diff_lte , by(xaxis)

save "LTE Example--Negative with Positive NPH.dta", replace


/*
*********************************************************************
*** Comparison of average hazard rate to true hazard
***
*** NOTE: do the following after you have run the Monte Carlo experiments!
***
*** Sims: 1,000
*** N: 1,000
*** Betas: (-3, 1)
*** Knot location: c(1, 4, 7)
*** X: uniformly drawn from [-2, 2]
*********************************************************************

foreach s of numlist 1(1)4 {
	local nv = S[`s',1]
	local pv = S[`s',2]
	local nm1v = S[`s',3]
	local nm2v = S[`s',4]
	local nwpv = S[`s',5]
	local pwnv = S[`s',6]

*** Negative
use "Negative_Hazard_`s'.dta", clear

local b = 1
foreach v of varlist TD_h CP_h S_h ASS_h {
	bys t: egen mnh`b' = mean(`v')
	local b = `b' + 1
}

duplicates drop t, force

preserve
	keep t true_h
	rename true_h mnh
	gen funcform = 5
	sort t funcform
	tempfile n_th 
	save `n_th', replace
restore

reshape long mnh, i(t) j(funcform)

append using `n_th'
gen scenario = 1

keep t mnh funcform scenario
sort t

tempfile n
save `n', replace


*** Positive
use "Positive_Hazard_`s'.dta", clear

local b = 1
foreach v of varlist TD_h CP_h S_h ASS_h {
	bys t: egen mnh`b' = mean(`v')
	local b = `b' + 1
}

duplicates drop t, force

preserve
	keep t true_h
	rename true_h mnh
	gen funcform = 5
	sort t funcform
	tempfile p_th 
	save `p_th', replace
restore

reshape long mnh, i(t) j(funcform)

append using `p_th'
gen scenario = 2

keep t mnh funcform scenario
sort t

tempfile p
save `p', replace


*** Non-Monotonic 1
use "NM1_Hazard_`s'.dta", clear

local b = 1
foreach v of varlist TD_h CP_h S_h ASS_h {
	bys t: egen mnh`b' = mean(`v')
	local b = `b' + 1
}

duplicates drop t, force

preserve
	keep t true_h
	rename true_h mnh
	gen funcform = 5
	sort t funcform
	tempfile nm1_th 
	save `nm1_th', replace
restore

reshape long mnh, i(t) j(funcform)

append using `nm1_th'
gen scenario = 3

keep t mnh funcform scenario
sort t

tempfile nm1
save `nm1', replace

*** Non-Monotonic 2
use "NM2_Hazard_`s'.dta", clear

local b = 1
foreach v of varlist TD_h CP_h S_h ASS_h {
	bys t: egen mnh`b' = mean(`v')
	local b = `b' + 1
}

duplicates drop t, force

preserve
	keep t true_h
	rename true_h mnh
	gen funcform = 5
	sort t funcform
	tempfile nm2_th 
	save `nm2_th', replace
restore

reshape long mnh, i(t) j(funcform)

append using `nm2_th'
gen scenario = 4

keep t mnh funcform scenario
sort t

tempfile nm2
save `nm2', replace

*** Negative with Positive NPH
use "NPH Positive_Hazard_`s'.dta", clear

local b = 1
foreach v of varlist TD_h CP_h S_h ASS_h {
	bys t: egen mnh`b' = mean(`v')
	local b = `b' + 1
}

duplicates drop t, force

preserve
	keep t true_h
	rename true_h mnh
	gen funcform = 5
	sort t funcform
	tempfile nwp_th 
	save `nwp_th', replace
restore

reshape long mnh, i(t) j(funcform)

append using `nwp_th'
gen scenario = 5

keep t mnh funcform scenario
sort t

tempfile nwp
save `nwp', replace


*** Positive with Negative NPH
use "NPH Negative_Hazard_`s'.dta", clear

local b = 1
foreach v of varlist TD_h CP_h S_h ASS_h {
	bys t: egen mnh`b' = mean(`v')
	local b = `b' + 1
}

duplicates drop t, force

preserve
	keep t true_h
	rename true_h mnh
	gen funcform = 5
	sort t funcform
	tempfile pwn_th 
	save `pwn_th', replace
restore

reshape long mnh, i(t) j(funcform)

append using `pwn_th'
gen scenario = 6

keep t mnh funcform scenario
sort t

tempfile pwn
save `pwn', replace

use `n', clear

append using `p'
append using `nm1'
append using `nm2'
append using `nwp'
append using `pwn'

sort scenario funcform t

lab def scenario_l 1 "Negative" 2 "Positive" 3 "Non-Monotonic 1" 4 "Non-Monotonic 2" 5 "Negative with Positive NPH" 6 "Positive with Negative NPH"
lab val scenario scenario_l

save "Average Hazard Rates--Scenario_`s'.dta", replace

}
*/


*************************************************************
*************************************************************
************************* Clare (2010) **********************
*************************************************************
*************************************************************

*************************************************************
*** Long-Term Effects Using Observed-Values Approach: Ideological Distance goes from its minimum to its maximum
*************************************************************
use "Clare.dta", clear

qui probit init_monad coalition dist_directional1 gov_ideology1 minority cap_1 prop_demborder1 prop_allyborder1 dep_trade veto_players peaceyears  peaceyears2 peaceyears3, robust
mat b = e(b)
mat V = e(V)

keep if e(sample) 
keep ccode1 year quarter init_monad coalition dist_directional1 gov_ideology1 minority cap_1 prop_demborder1 prop_allyborder1 dep_trade veto_players peaceyears  peaceyears2 peaceyears3

tempfile tdata
save `tdata', replace

set seed 648
drawnorm b_coalition b_dist_directional1 b_gov_ideology1 b_minority b_cap_1 b_prop_demborder1 b_prop_allyborder1 b_dep_trade b_veto_players b_peaceyears b_peaceyears2 b_peaceyears3 b_cons, /*
*/	mean(b) cov(V) n(1000) clear

merge using `tdata'

* Rows: x (permanent), t0, t1
* Columns: Baseline, counterfactual, then simulations
mat S = (-110, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120,120,120,120,120,120,120,120,120,120,120,120,120,120,120,120 \ /*
*/		 5, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30 \ /*
*/		 5, 5, 0, 1, 2, 3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24 )

mat SIMS = S[.,3...]
local simno = colsof(SIMS)
local totalsim = colsof(S)

preserve
	replace dist_directional1 = -110
	predict p_min
	replace dist_directional1 = 120
	predict p_max
	
	rename peaceyears t
	keep t p_min p_max
	save "ISPP.dta", replace
restore

tempname alte
postfile `alte' xaxis time0 time1 p0 p0_lo p0_hi p1 p1_lo p1_hi lte lte_lo lte_hi /*
*/	using "Clare--LTE--Average.dta", replace

qui foreach c of numlist 1(1)`simno' {
	local x = S[1,`c']
	local t0 = S[2,`c']
	local t1 = S[3,`c']

	tempvar pr0_`c' pr1_`c' lte_`c'
	gen `pr0_`c'' = .
	gen `pr1_`c'' = .
	gen `lte_`c'' = .

	qui forvalues s = 1/1000 {
		tempvar lte_s`s'
		foreach v of numlist 0 1 {
			tempvar pr`v'_s`s'
			gen `pr`v'_s`s'' = normal(b_cons[`s'] + coalition*b_coalition[`s'] + `x'*b_dist_directional1[`s'] + gov_ideology1*b_gov_ideology1[`s'] +  minority*b_minority[`s'] + cap_1*b_cap_1[`s'] + prop_demborder1*b_prop_demborder1[`s'] +  prop_allyborder1*b_prop_allyborder1[`s'] + dep_trade*b_dep_trade[`s'] + veto_players*b_veto_players[`s'] + `t`v''*b_peaceyears[`s'] + `t`v''^2*b_peaceyears2[`s'] + `t`v''^3*b_peaceyears3[`s'])
			sum `pr`v'_s`s'', meanonly
			replace `pr`v'_`c'' = r(mean) in `s'
		}
		gen `lte_s`s'' = `pr1_s`s'' - `pr0_s`s''
		sum `lte_s`s'', meanonly
		replace `lte_`c'' = r(mean) in `s'
	
		if mod(`s',100) == 0 {
			nois display "." _c
			if mod(`s',1000) == 0 {
				nois display "Simulation #`c' of `totalsim'"
			}
		}	
	}
	sum `pr0_`c'', meanonly
	local p0 = r(mean)
	_pctile `pr0_`c'', perc(2.5 97.5)
	local p0_lo = r(r1)
	local p0_hi = r(r2)
	
	sum `pr1_`c'', meanonly
	local p1 = r(mean)
	_pctile `pr1_`c'', perc(2.5 97.5)
	local p1_lo = r(r1)
	local p1_hi = r(r2)

	sum `lte_`c'', meanonly
	local lte = r(mean)
	_pctile `lte_`c'', perc(2.5 97.5)
	local lte_lo = r(r1)
	local lte_hi = r(r2)

	local xaxis = `c' - 2
	
	cap drop __*
	post `alte' (`xaxis') (`t0') (`t1') (`p0') (`p0_lo') (`p0_hi') (`p1') (`p1_lo') (`p1_hi') (`lte') (`lte_lo') (`lte_hi')	
}

postclose `alte'


*************************************************************
*** Long-Term Effects for Four Simulation Scenarios as Ideological Distance Increases from its Minimum to its Maximum
*************************************************************
use "Clare.dta", clear

qui probit init_monad coalition dist_directional1 gov_ideology1 minority cap_1 prop_demborder1 prop_allyborder1 dep_trade veto_players peaceyears  peaceyears2 peaceyears3, robust
keep if e(sample) 

set seed 648
estsimp probit init_monad coalition dist_directional1 gov_ideology1 minority cap_1 prop_demborder1 prop_allyborder1 dep_trade veto_players peaceyears  peaceyears2 peaceyears3, robust

global no "coalition 0 dist_directional1 min gov_ideology1 min minority 1 cap_1 p25 prop_demborder1 p75 prop_allyborder1 p75 dep_trade p75 veto_players p75 peaceyears 5 peaceyears2 25 peaceyears3 125"
global low "coalition 1 dist_directional1 p75 gov_ideology1 p75 minority 0 cap_1 p75 prop_demborder1 mean prop_allyborder1 mean dep_trade min veto_players p25 peaceyears 5 peaceyears2 25 peaceyears3 125"
global mod "coalition 1 dist_directional1 p95 gov_ideology1 p80 minority 0 cap_1 p90 prop_demborder1 p5 prop_allyborder1 p5 dep_trade min veto_players p5 peaceyears 5 peaceyears2 25 peaceyears3 125"
global hi "coalition 1 dist_directional1 200 gov_ideology1 200 minority 0 cap_1 max prop_demborder1 min prop_allyborder1 min dep_trade min veto_players min peaceyears 5 peaceyears2 25 peaceyears3 125"

* Rows: x (1-time), x (temporary), x (permanent), t0_0, t0_1, t1_0, t1_1
* Columns: Baseline, counterfactual, then simulations
mat S = (-110, 120,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0 \ /*
*/  	 -110, 120, 120, 120, 120,  0,  0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0 \ /*
*/		 -110, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120, 120,120,120,120,120,120 \ /*
*/		 5, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20 \ /*
*/		 5, 5, 0, 1, 2, 3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14 )

mat SIMS = S[.,3...]
local simno = colsof(SIMS)

* No risk
setx $no
simqi, prval(1)

* Low risk:
setx $low
simqi, prval(1)

* Moderate risk: 0.10
setx $mod
simqi, prval(1)

* High (coin flip) risk:
setx $hi
simqi, prval(1)

tempname slte
postfile `slte' xaxis time0 time1 str20 scenario /*
*/	x1 pr_s1_0 pr_s1_0_lo pr_s1_0_hi pr_s1_1 pr_s1_1_lo pr_s1_1_hi lte_s1 lte_s1_lo lte_s1_hi /*
*/	x2 pr_s2_0 pr_s2_0_lo pr_s2_0_hi pr_s2_1 pr_s2_1_lo pr_s2_1_hi lte_s2 lte_s2_lo lte_s2_hi /*
*/	x3 pr_s3_0 pr_s3_0_lo pr_s3_0_hi pr_s3_1 pr_s3_1_lo pr_s3_1_hi lte_s3 lte_s3_lo lte_s3_hi /*
*/	using "Clare--LTE--Scenarios.dta", replace

foreach c of numlist 1(1)`simno' {
	local x1 = S[1,`c']
	local x2 = S[2,`c']
	local x3 = S[3,`c']
	local t0 = S[4,`c']
	local t1 = S[5,`c']
	
	qui foreach i in no low mod hi {
		foreach s of numlist 1(1)3 {
			foreach o of numlist 0 1 {
				tempvar prob_s`s'_`o'
				setx $`i'
				setx peaceyears `t`o'' peaceyears2 `t`o''^2 peaceyears3 `t`o''^3 dist_directional1 `x`s''
				simqi, prval(1) genpr(`prob_s`s'_`o'')
				sum `prob_s`s'_`o'', meanonly
				local pr_s`s'_`o' = r(mean)
				_pctile `prob_s`s'_`o'', perc(2.5 97.5)
				local pr_s`s'_`o'_lo = r(r1)
				local pr_s`s'_`o'_hi = r(r2)
			}	
	
			tempvar lteff_s`s'
			gen `lteff_s`s'' = `prob_s`s'_1' - `prob_s`s'_0'
			sum `lteff_s`s'', meanonly
			local lte_s`s' = r(mean)
			_pctile `lteff_s`s'', perc(2.5 97.5)
			local lte_s`s'_lo = r(r1)
			local lte_s`s'_hi = r(r2)
		}
		
		local xaxis = `c' - 2
	
		post `slte' (`xaxis') (`t0') (`t1') ("`i'") /*
		*/	(`x1') (`pr_s1_0') (`pr_s1_0_lo') (`pr_s1_0_hi') (`pr_s1_1') (`pr_s1_1_lo') (`pr_s1_1_hi') (`lte_s1') (`lte_s1_lo') (`lte_s1_hi') /*
		*/	(`x2') (`pr_s2_0') (`pr_s2_0_lo') (`pr_s2_0_hi') (`pr_s2_1') (`pr_s2_1_lo') (`pr_s2_1_hi') (`lte_s2') (`lte_s2_lo') (`lte_s2_hi') /*
		*/	(`x3') (`pr_s3_0') (`pr_s3_0_lo') (`pr_s3_0_hi') (`pr_s3_1') (`pr_s3_1_lo') (`pr_s3_1_hi') (`lte_s3') (`lte_s3_lo') (`lte_s3_hi')			
	}	
	
	if mod(`c',1) == 0 {
		nois display "." _c
	}
}

postclose `slte'



*****************************************************************
*****************************************************************
*********************** APPENDIX ********************************
*****************************************************************
*****************************************************************

*****************************************************************
*** Four functional forms of duration dependence
***
*** Original data sets are created in R ("Generate Functional Form Examples.R")
*****************************************************************
*** Negative
use "NegativeFF.dta", clear

btscs y tim group, g(t) nspline(3)
replace t = (t + 1)
gen t2 = t^2
gen t3 = t^3

set seed 648
estsimp logit y x t t2 t3

qui foreach k in h_t h h_lo h_hi {
	gen `k' = .
}

local b = 1
setx x 0
qui foreach v of numlist 1(1)25 {
	setx t `v' t2 `v'^2 t3 `v'^3
	simqi, prval(1) genpr(prh)
	sum prh, meanonly
	local h = r(mean)
	_pctile prh, perc(2.5 97.5)
	local h_lo = r(r1)
	local h_hi = r(r2)
	
	replace h_t = `v' in `b'
	replace h = `h' in `b'
	replace h_lo = `h_lo' in `b'
	replace h_hi = `h_hi' in `b'
	
	drop prh
	local b = `b' + 1
}

gen scen = "Negative"

tempfile neg
save `neg', replace

*** Positive
use "PositiveFF.dta", clear

btscs y tim group, g(t) nspline(3)
replace t = (t + 1)
gen t2 = t^2
gen t3 = t^3

set seed 648
estsimp logit y x t t2 t3

qui foreach k in h_t h h_lo h_hi {
	gen `k' = .
}

local b = 1
setx x 0
qui foreach v of numlist 1(1)25 {
	setx t `v' t2 `v'^2 t3 `v'^3
	simqi, prval(1) genpr(prh)
	sum prh, meanonly
	local h = r(mean)
	_pctile prh, perc(2.5 97.5)
	local h_lo = r(r1)
	local h_hi = r(r2)
	
	replace h_t = `v' in `b'
	replace h = `h' in `b'
	replace h_lo = `h_lo' in `b'
	replace h_hi = `h_hi' in `b'
	
	drop prh
	local b = `b' + 1
}

gen scen = "Positive"

tempfile pos
save `pos', replace

*** Non-Monotonic duration dependence 1
use "Non-Monotonic1FF.dta", clear

btscs y tim group, g(t) nspline(3)
replace t = (t + 1)
gen t2 = t^2
gen t3 = t^3

set seed 648
estsimp logit y x t t2 t3

qui foreach k in h_t h h_lo h_hi {
	gen `k' = .
}

local b = 1
setx x 0
qui foreach v of numlist 1(1)25 {
	setx t `v' t2 `v'^2 t3 `v'^3
	simqi, prval(1) genpr(prh)
	sum prh, meanonly
	local h = r(mean)
	_pctile prh, perc(2.5 97.5)
	local h_lo = r(r1)
	local h_hi = r(r2)
	
	replace h_t = `v' in `b'
	replace h = `h' in `b'
	replace h_lo = `h_lo' in `b'
	replace h_hi = `h_hi' in `b'
	
	drop prh
	local b = `b' + 1
}

gen scen = "Non-Monotonic 1"

tempfile nm1
save `nm1', replace


*** Non-Monotonic duration dependence 2
use "Non-Monotonic2FF.dta", clear

btscs y tim group, g(t) nspline(3)
replace t = (t + 1)
gen t2 = t^2
gen t3 = t^3

set seed 648
estsimp logit y x t t2 t3

qui foreach k in h_t h h_lo h_hi {
	gen `k' = .
}

local b = 1
setx x 0
qui foreach v of numlist 1(1)25 {
	setx t `v' t2 `v'^2 t3 `v'^3
	simqi, prval(1) genpr(prh)
	sum prh, meanonly
	local h = r(mean)
	_pctile prh, perc(2.5 97.5)
	local h_lo = r(r1)
	local h_hi = r(r2)
	
	replace h_t = `v' in `b'
	replace h = `h' in `b'
	replace h_lo = `h_lo' in `b'
	replace h_hi = `h_hi' in `b'
	
	drop prh
	local b = `b' + 1
}

gen scen = "Non-Monotonic 2"

append using `neg'
append using `pos'
append using `nm1'

keep h_t h h_lo h_hi scen
drop if h_t == .

save "Functional Forms.dta", replace


*****************************************************************
*** Long-term effects for four scenarios for negative duration dependence
*****************************************************************
use "Negative.dta", clear

btscs y tim group, g(t) nspline(3)
replace t = (t + 1)
gen t2 = t^2
gen t3 = t^3

set seed 648
estsimp logit y x t t2 t3

************************ All Types of Effects *****************************
* Rows: x (1-time), t0_1, t0_2, t0_3, t0_4, t1
* Columns: Baseline, counterfactual, then simulations
mat S = (0 , .5,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 \ /*
*/		 0,   0,  1, 2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25 \ /*
*/		 4,   4,  5, 6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29 \ /*
*/		 9,   9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34 \ /*
*/		 19, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44 \ /*
*/		 5,  5,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9,  10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24)

mat SIMS = S[.,3...]
local simno = colsof(SIMS)

tempname lteneg
postfile `lteneg' xaxis time0_1 time0_2 time0_3 time0_4 time1 /*
*/	x1 pr_s1_0 pr_s1_0_lo pr_s1_0_hi pr_s1_1 pr_s1_1_lo pr_s1_1_hi lte_s1 lte_s1_lo lte_s1_hi /*
*/	pr_s2_0 pr_s2_0_lo pr_s2_0_hi pr_s2_1 pr_s2_1_lo pr_s2_1_hi lte_s2 lte_s2_lo lte_s2_hi /*
*/	pr_s3_0 pr_s3_0_lo pr_s3_0_hi pr_s3_1 pr_s3_1_lo pr_s3_1_hi lte_s3 lte_s3_lo lte_s3_hi /*
*/	pr_s4_0 pr_s4_0_lo pr_s4_0_hi pr_s4_1 pr_s4_1_lo pr_s4_1_hi lte_s4 lte_s4_lo lte_s4_hi /*
*/	using "LTE Example--Scenarios--Negative.dta", replace


qui foreach c of numlist 1(1)`simno' {
	local x1 = S[1,`c']
	local t0_1 = S[2,`c']
	local t0_2 = S[3,`c']
	local t0_3 = S[4,`c']
	local t0_4 = S[5,`c']
	local t1 = S[6,`c']

	foreach s of numlist 1(1)4 {
		tempvar prob 
		setx x `x1' t `t1' t2 `t1'^2 t3 `t1'^3
		simqi, prval(1) genpr(`prob')
		sum `prob', meanonly
		local pr_s`s'_1 = r(mean)
		_pctile `prob', perc(2.5 97.5)
		local pr_s`s'_1_lo = r(r1)
		local pr_s`s'_1_hi = r(r2)
		
		tempvar prob_s`s'_0
		setx x `x1' t `t0_`s'' t2 `t0_`s''^2 t3 `t0_`s''^3
				
		simqi, prval(1) genpr(`prob_s`s'_0')
		sum `prob_s`s'_0', meanonly
		local pr_s`s'_0 = r(mean)
		_pctile `prob_s`s'_0', perc(2.5 97.5)
		local pr_s`s'_0_lo = r(r1)
		local pr_s`s'_0_hi = r(r2)
	
		tempvar lteff_s`s'
		gen `lteff_s`s'' = `prob' - `prob_s`s'_0'
		sum `lteff_s`s'', meanonly
		local lte_s`s' = r(mean)
		_pctile `lteff_s`s'', perc(2.5 97.5)
		local lte_s`s'_lo = r(r1)
		local lte_s`s'_hi = r(r2)
	}
	
	local xaxis = `c' - 2
	
	post `lteneg' (`xaxis') (`t0_1') (`t0_2') (`t0_3') (`t0_4') (`t1') /*
*/	(`x1') (`pr_s1_0') (`pr_s1_0_lo') (`pr_s1_0_hi') (`pr_s1_1') (`pr_s1_1_lo') (`pr_s1_1_hi') (`lte_s1') (`lte_s1_lo') (`lte_s1_hi') /*
*/	(`pr_s2_0') (`pr_s2_0_lo') (`pr_s2_0_hi') (`pr_s2_1') (`pr_s2_1_lo') (`pr_s2_1_hi') (`lte_s2') (`lte_s2_lo') (`lte_s2_hi') /*
*/	(`pr_s3_0') (`pr_s3_0_lo') (`pr_s3_0_hi') (`pr_s3_1') (`pr_s3_1_lo') (`pr_s3_1_hi') (`lte_s3') (`lte_s3_lo') (`lte_s3_hi') /*
*/	(`pr_s4_0') (`pr_s4_0_lo') (`pr_s4_0_hi') (`pr_s4_1') (`pr_s4_1_lo') (`pr_s4_1_hi') (`lte_s4') (`lte_s4_lo') (`lte_s4_hi') 
}

postclose `lteneg'


*****************************************************************
*** Long-term effects for four scenarios for positive duration dependence
*****************************************************************
use "Positive.dta", clear

btscs y tim group, g(t) nspline(3)
replace t = (t + 1)
gen t2 = t^2
gen t3 = t^3

set seed 648
estsimp logit y x t t2 t3

************************ All Types of Effects *****************************
* Rows: x (1-time), t0_1, t0_2, t0_3, t0_4, t1
* Columns: Baseline, counterfactual, then simulations
mat S = (0 , .5,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 \ /*
*/		 0,   0,  1, 2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25 \ /*
*/		 4,   4,  5, 6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29 \ /*
*/		 9,   9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34 \ /*
*/		 19, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44 \ /*
*/		 5,  5,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9,  10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24)

mat SIMS = S[.,3...]
local simno = colsof(SIMS)

tempname ltepos
postfile `ltepos' xaxis time0_1 time0_2 time0_3 time0_4 time1 /*
*/	x1 pr_s1_0 pr_s1_0_lo pr_s1_0_hi pr_s1_1 pr_s1_1_lo pr_s1_1_hi lte_s1 lte_s1_lo lte_s1_hi /*
*/	pr_s2_0 pr_s2_0_lo pr_s2_0_hi pr_s2_1 pr_s2_1_lo pr_s2_1_hi lte_s2 lte_s2_lo lte_s2_hi /*
*/	pr_s3_0 pr_s3_0_lo pr_s3_0_hi pr_s3_1 pr_s3_1_lo pr_s3_1_hi lte_s3 lte_s3_lo lte_s3_hi /*
*/	pr_s4_0 pr_s4_0_lo pr_s4_0_hi pr_s4_1 pr_s4_1_lo pr_s4_1_hi lte_s4 lte_s4_lo lte_s4_hi /*
*/	using "LTE Example--Scenarios--Positive.dta", replace


qui foreach c of numlist 1(1)`simno' {
	local x1 = S[1,`c']
	local t0_1 = S[2,`c']
	local t0_2 = S[3,`c']
	local t0_3 = S[4,`c']
	local t0_4 = S[5,`c']
	local t1 = S[6,`c']

	foreach s of numlist 1(1)4 {
		tempvar prob 
		setx x `x1' t `t1' t2 `t1'^2 t3 `t1'^3
		simqi, prval(1) genpr(`prob')
		sum `prob', meanonly
		local pr_s`s'_1 = r(mean)
		_pctile `prob', perc(2.5 97.5)
		local pr_s`s'_1_lo = r(r1)
		local pr_s`s'_1_hi = r(r2)
		
		tempvar prob_s`s'_0
		setx x `x1' t `t0_`s'' t2 `t0_`s''^2 t3 `t0_`s''^3
				
		simqi, prval(1) genpr(`prob_s`s'_0')
		sum `prob_s`s'_0', meanonly
		local pr_s`s'_0 = r(mean)
		_pctile `prob_s`s'_0', perc(2.5 97.5)
		local pr_s`s'_0_lo = r(r1)
		local pr_s`s'_0_hi = r(r2)
	
		tempvar lteff_s`s'
		gen `lteff_s`s'' = `prob' - `prob_s`s'_0'
		sum `lteff_s`s'', meanonly
		local lte_s`s' = r(mean)
		_pctile `lteff_s`s'', perc(2.5 97.5)
		local lte_s`s'_lo = r(r1)
		local lte_s`s'_hi = r(r2)
	}
	
	local xaxis = `c' - 2
	
	post `ltepos' (`xaxis') (`t0_1') (`t0_2') (`t0_3') (`t0_4') (`t1') /*
*/	(`x1') (`pr_s1_0') (`pr_s1_0_lo') (`pr_s1_0_hi') (`pr_s1_1') (`pr_s1_1_lo') (`pr_s1_1_hi') (`lte_s1') (`lte_s1_lo') (`lte_s1_hi') /*
*/	(`pr_s2_0') (`pr_s2_0_lo') (`pr_s2_0_hi') (`pr_s2_1') (`pr_s2_1_lo') (`pr_s2_1_hi') (`lte_s2') (`lte_s2_lo') (`lte_s2_hi') /*
*/	(`pr_s3_0') (`pr_s3_0_lo') (`pr_s3_0_hi') (`pr_s3_1') (`pr_s3_1_lo') (`pr_s3_1_hi') (`lte_s3') (`lte_s3_lo') (`lte_s3_hi') /*
*/	(`pr_s4_0') (`pr_s4_0_lo') (`pr_s4_0_hi') (`pr_s4_1') (`pr_s4_1_lo') (`pr_s4_1_hi') (`lte_s4') (`lte_s4_lo') (`lte_s4_hi') 
}

postclose `ltepos'


*****************************************************************
*** Long-term effects for four scenarios for non-monotonic 1 duration dependence
*****************************************************************
use "NM1.dta", clear

btscs y tim group, g(t) nspline(3)
replace t = (t + 1)
gen t2 = t^2
gen t3 = t^3

set seed 648
estsimp logit y x t t2 t3

************************ All Types of Effects *****************************
* Rows: x (1-time), t0_1, t0_2, t0_3, t0_4, t1
* Columns: Baseline, counterfactual, then simulations
mat S = (0 , .5,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 \ /*
*/		 0,   0,  1, 2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25 \ /*
*/		 4,   4,  5, 6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29 \ /*
*/		 9,   9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34 \ /*
*/		 19, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44 \ /*
*/		 5,  5,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9,  10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24)

mat SIMS = S[.,3...]
local simno = colsof(SIMS)

tempname ltenm1
postfile `ltenm1' xaxis time0_1 time0_2 time0_3 time0_4 time1 /*
*/	x1 pr_s1_0 pr_s1_0_lo pr_s1_0_hi pr_s1_1 pr_s1_1_lo pr_s1_1_hi lte_s1 lte_s1_lo lte_s1_hi /*
*/	pr_s2_0 pr_s2_0_lo pr_s2_0_hi pr_s2_1 pr_s2_1_lo pr_s2_1_hi lte_s2 lte_s2_lo lte_s2_hi /*
*/	pr_s3_0 pr_s3_0_lo pr_s3_0_hi pr_s3_1 pr_s3_1_lo pr_s3_1_hi lte_s3 lte_s3_lo lte_s3_hi /*
*/	pr_s4_0 pr_s4_0_lo pr_s4_0_hi pr_s4_1 pr_s4_1_lo pr_s4_1_hi lte_s4 lte_s4_lo lte_s4_hi /*
*/	using "LTE Example--Scenarios--NM1.dta", replace


qui foreach c of numlist 1(1)`simno' {
	local x1 = S[1,`c']
	local t0_1 = S[2,`c']
	local t0_2 = S[3,`c']
	local t0_3 = S[4,`c']
	local t0_4 = S[5,`c']
	local t1 = S[6,`c']

	foreach s of numlist 1(1)4 {
		tempvar prob 
		setx x `x1' t `t1' t2 `t1'^2 t3 `t1'^3
		simqi, prval(1) genpr(`prob')
		sum `prob', meanonly
		local pr_s`s'_1 = r(mean)
		_pctile `prob', perc(2.5 97.5)
		local pr_s`s'_1_lo = r(r1)
		local pr_s`s'_1_hi = r(r2)
		
		tempvar prob_s`s'_0
		setx x `x1' t `t0_`s'' t2 `t0_`s''^2 t3 `t0_`s''^3
				
		simqi, prval(1) genpr(`prob_s`s'_0')
		sum `prob_s`s'_0', meanonly
		local pr_s`s'_0 = r(mean)
		_pctile `prob_s`s'_0', perc(2.5 97.5)
		local pr_s`s'_0_lo = r(r1)
		local pr_s`s'_0_hi = r(r2)
	
		tempvar lteff_s`s'
		gen `lteff_s`s'' = `prob' - `prob_s`s'_0'
		sum `lteff_s`s'', meanonly
		local lte_s`s' = r(mean)
		_pctile `lteff_s`s'', perc(2.5 97.5)
		local lte_s`s'_lo = r(r1)
		local lte_s`s'_hi = r(r2)
	}
	
	local xaxis = `c' - 2
	
	post `ltenm1' (`xaxis') (`t0_1') (`t0_2') (`t0_3') (`t0_4') (`t1') /*
*/	(`x1') (`pr_s1_0') (`pr_s1_0_lo') (`pr_s1_0_hi') (`pr_s1_1') (`pr_s1_1_lo') (`pr_s1_1_hi') (`lte_s1') (`lte_s1_lo') (`lte_s1_hi') /*
*/	(`pr_s2_0') (`pr_s2_0_lo') (`pr_s2_0_hi') (`pr_s2_1') (`pr_s2_1_lo') (`pr_s2_1_hi') (`lte_s2') (`lte_s2_lo') (`lte_s2_hi') /*
*/	(`pr_s3_0') (`pr_s3_0_lo') (`pr_s3_0_hi') (`pr_s3_1') (`pr_s3_1_lo') (`pr_s3_1_hi') (`lte_s3') (`lte_s3_lo') (`lte_s3_hi') /*
*/	(`pr_s4_0') (`pr_s4_0_lo') (`pr_s4_0_hi') (`pr_s4_1') (`pr_s4_1_lo') (`pr_s4_1_hi') (`lte_s4') (`lte_s4_lo') (`lte_s4_hi') 
}

postclose `ltenm1'


*****************************************************************
*** Long-term effects for four scenarios for non-monotonic 2 duration dependence
*****************************************************************
use "NM2.dta", clear

btscs y tim group, g(t) nspline(3)
replace t = (t + 1)
gen t2 = t^2
gen t3 = t^3

set seed 648
estsimp logit y x t t2 t3

************************ All Types of Effects *****************************
* Rows: x (1-time), t0_1, t0_2, t0_3, t0_4, t1
* Columns: Baseline, counterfactual, then simulations
mat S = (0 , .5,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 \ /*
*/		 0,   0,  1, 2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25 \ /*
*/		 4,   4,  5, 6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29 \ /*
*/		 9,   9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34 \ /*
*/		 19, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44 \ /*
*/		 5,  5,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9,  10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24)

mat SIMS = S[.,3...]
local simno = colsof(SIMS)

tempname ltenm2
postfile `ltenm2' xaxis time0_1 time0_2 time0_3 time0_4 time1 /*
*/	x1 pr_s1_0 pr_s1_0_lo pr_s1_0_hi pr_s1_1 pr_s1_1_lo pr_s1_1_hi lte_s1 lte_s1_lo lte_s1_hi /*
*/	pr_s2_0 pr_s2_0_lo pr_s2_0_hi pr_s2_1 pr_s2_1_lo pr_s2_1_hi lte_s2 lte_s2_lo lte_s2_hi /*
*/	pr_s3_0 pr_s3_0_lo pr_s3_0_hi pr_s3_1 pr_s3_1_lo pr_s3_1_hi lte_s3 lte_s3_lo lte_s3_hi /*
*/	pr_s4_0 pr_s4_0_lo pr_s4_0_hi pr_s4_1 pr_s4_1_lo pr_s4_1_hi lte_s4 lte_s4_lo lte_s4_hi /*
*/	using "LTE Example--Scenarios--NM2.dta", replace


qui foreach c of numlist 1(1)`simno' {
	local x1 = S[1,`c']
	local t0_1 = S[2,`c']
	local t0_2 = S[3,`c']
	local t0_3 = S[4,`c']
	local t0_4 = S[5,`c']
	local t1 = S[6,`c']

	foreach s of numlist 1(1)4 {
		tempvar prob 
		setx x `x1' t `t1' t2 `t1'^2 t3 `t1'^3
		simqi, prval(1) genpr(`prob')
		sum `prob', meanonly
		local pr_s`s'_1 = r(mean)
		_pctile `prob', perc(2.5 97.5)
		local pr_s`s'_1_lo = r(r1)
		local pr_s`s'_1_hi = r(r2)
		
		tempvar prob_s`s'_0
		setx x `x1' t `t0_`s'' t2 `t0_`s''^2 t3 `t0_`s''^3
				
		simqi, prval(1) genpr(`prob_s`s'_0')
		sum `prob_s`s'_0', meanonly
		local pr_s`s'_0 = r(mean)
		_pctile `prob_s`s'_0', perc(2.5 97.5)
		local pr_s`s'_0_lo = r(r1)
		local pr_s`s'_0_hi = r(r2)
	
		tempvar lteff_s`s'
		gen `lteff_s`s'' = `prob' - `prob_s`s'_0'
		sum `lteff_s`s'', meanonly
		local lte_s`s' = r(mean)
		_pctile `lteff_s`s'', perc(2.5 97.5)
		local lte_s`s'_lo = r(r1)
		local lte_s`s'_hi = r(r2)
	}
	
	local xaxis = `c' - 2
	
	post `ltenm2' (`xaxis') (`t0_1') (`t0_2') (`t0_3') (`t0_4') (`t1') /*
*/	(`x1') (`pr_s1_0') (`pr_s1_0_lo') (`pr_s1_0_hi') (`pr_s1_1') (`pr_s1_1_lo') (`pr_s1_1_hi') (`lte_s1') (`lte_s1_lo') (`lte_s1_hi') /*
*/	(`pr_s2_0') (`pr_s2_0_lo') (`pr_s2_0_hi') (`pr_s2_1') (`pr_s2_1_lo') (`pr_s2_1_hi') (`lte_s2') (`lte_s2_lo') (`lte_s2_hi') /*
*/	(`pr_s3_0') (`pr_s3_0_lo') (`pr_s3_0_hi') (`pr_s3_1') (`pr_s3_1_lo') (`pr_s3_1_hi') (`lte_s3') (`lte_s3_lo') (`lte_s3_hi') /*
*/	(`pr_s4_0') (`pr_s4_0_lo') (`pr_s4_0_hi') (`pr_s4_1') (`pr_s4_1_lo') (`pr_s4_1_hi') (`lte_s4') (`lte_s4_lo') (`lte_s4_hi') 
}

postclose `ltenm2'


*****************************************************************
*** Variations of long-term effects based on the counterfactual of interest for Positive duration dependence
*****************************************************************
use "Positive.dta", clear

btscs y tim group, g(t) nspline(3)
replace t = (t + 1)
gen t2 = t^2
gen t3 = t^3

set seed 648
estsimp logit y x t t2 t3

************************ All Types of Effects *****************************
* Rows: x (1-time), x (temporary), x (permanent), t0_0, t0_1, t1_0, t1_1
* Columns: Baseline, counterfactual, then simulations
mat S = (0 , .5,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0 \ /*
*/  	 0 , .5, .5, .5, .5,  0,  0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0 \ /*
*/		 0 , .5, .5, .5, .5, .5, .5, .5, .5, .5, .5, .5,.5,.5,.5,.5,.5 \ /*
*/		 13, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28 \ /*
*/		 13, 13,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14 \ /*
*/		 13, 13,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14 \ /*
*/		 13, 13,  0,  1,  2,  3,  4,  0,  1,  2,  3,  4,  5, 6, 7, 8, 9)

mat SIMS = S[.,3...]
local simno = colsof(SIMS)


tempname lteneg
postfile `lteneg' xaxis time0_0 time0_1 time1_0 time1_1 /*
*/	x1 pr_s1_0_0 pr_s1_0_0_lo pr_s1_0_0_hi pr_s1_0_1 pr_s1_0_1_lo pr_s1_0_1_hi pr_s1_1_0 pr_s1_1_0_lo pr_s1_1_0_hi pr_s1_1_1 pr_s1_1_1_lo pr_s1_1_1_hi lte_s1 lte_s1_lo lte_s1_hi /*
*/	x2 pr_s2_0_0 pr_s2_0_0_lo pr_s2_0_0_hi pr_s2_0_1 pr_s2_0_1_lo pr_s2_0_1_hi pr_s2_1_0 pr_s2_1_0_lo pr_s2_1_0_hi pr_s2_1_1 pr_s2_1_1_lo pr_s2_1_1_hi lte_s2 lte_s2_lo lte_s2_hi /*
*/	x3 pr_s3_0_0 pr_s3_0_0_lo pr_s3_0_0_hi pr_s3_0_1 pr_s3_0_1_lo pr_s3_0_1_hi pr_s3_1_0 pr_s3_1_0_lo pr_s3_1_0_hi pr_s3_1_1 pr_s3_1_1_lo pr_s3_1_1_hi lte_s3 lte_s3_lo lte_s3_hi /*
*/	using "LTE Example--Positive.dta", replace


qui foreach c of numlist 1(1)`simno' {
	local x1 = S[1,`c']
	local x2 = S[2,`c']
	local x3 = S[3,`c']
	local t0_0 = S[4,`c']
	local t0_1 = S[5,`c']
	local t1_0 = S[6,`c']
	local t1_1 = S[7,`c']
	
	*** Y_{t-1} == 0 & Y_{t-1} == 1	
	foreach s of numlist 1(1)3 {
		foreach b of numlist 0 1 {
			foreach o of numlist 0 1 {
				tempvar prob_s`s'_`o'_`b'
				setx x `x`s'' t `t`o'_`b'' t2 `t`o'_`b''^2 t3 `t`o'_`b''^3
				simqi, prval(1) genpr(`prob_s`s'_`o'_`b'')
				sum `prob_s`s'_`o'_`b'', meanonly
				local pr_s`s'_`o'_`b' = r(mean)
				_pctile `prob_s`s'_`o'_`b'', perc(2.5 97.5)
				local pr_s`s'_`o'_`b'_lo = r(r1)
				local pr_s`s'_`o'_`b'_hi = r(r2)
			}
		}
	
		tempvar lteff_s`s'
		gen `lteff_s`s'' = `prob_s`s'_0_1' - `prob_s`s'_0_0'
		sum `lteff_s`s'', meanonly
		local lte_s`s' = r(mean)
		_pctile `lteff_s`s'', perc(2.5 97.5)
		local lte_s`s'_lo = r(r1)
		local lte_s`s'_hi = r(r2)
	}
	
	local xaxis = `c' - 2
	
	post `lteneg' (`xaxis') (`t0_0') (`t0_1') (`t1_0') (`t1_1') /*
*/	(`x1') (`pr_s1_0_0') (`pr_s1_0_0_lo') (`pr_s1_0_0_hi') (`pr_s1_0_1') (`pr_s1_0_1_lo') (`pr_s1_0_1_hi') (`pr_s1_1_0') (`pr_s1_1_0_lo') (`pr_s1_1_0_hi') (`pr_s1_1_1') (`pr_s1_1_1_lo') (`pr_s1_1_1_hi') (`lte_s1') (`lte_s1_lo') (`lte_s1_hi') /*
*/	(`x2') (`pr_s2_0_0') (`pr_s2_0_0_lo') (`pr_s2_0_0_hi') (`pr_s2_0_1') (`pr_s2_0_1_lo') (`pr_s2_0_1_hi') (`pr_s2_1_0') (`pr_s2_1_0_lo') (`pr_s2_1_0_hi') (`pr_s2_1_1') (`pr_s2_1_1_lo') (`pr_s2_1_1_hi') (`lte_s2') (`lte_s2_lo') (`lte_s2_hi') /*
*/	(`x3') (`pr_s3_0_0') (`pr_s3_0_0_lo') (`pr_s3_0_0_hi') (`pr_s3_0_1') (`pr_s3_0_1_lo') (`pr_s3_0_1_hi') (`pr_s3_1_0') (`pr_s3_1_0_lo') (`pr_s3_1_0_hi') (`pr_s3_1_1') (`pr_s3_1_1_lo') (`pr_s3_1_1_hi') (`lte_s3') (`lte_s3_lo') (`lte_s3_hi')
}

postclose `lteneg'




*****************************************************************
*** Variations of long-term effects based on the counterfactual of interest for non-monotonic 1 duration dependence
*****************************************************************
use "NM1.dta", clear

btscs y tim group, g(t) nspline(3)
replace t = (t + 1)
gen t2 = t^2
gen t3 = t^3

set seed 648
*rename gg x
estsimp logit y x t t2 t3

************************ All Types of Effects *****************************
* Rows: x (1-time), x (temporary), x (permanent), t0_0, t0_1, t1_0, t1_1
* Columns: Baseline, counterfactual, then simulations
mat S = (0 , .5,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0 \ /*
*/  	 0 , .5, .5, .5, .5,  0,  0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0 \ /*
*/		 0 , .5, .5, .5, .5, .5, .5, .5, .5, .5, .5, .5,.5,.5,.5,.5,.5 \ /*
*/		 13, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28 \ /*
*/		 13, 13,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14 \ /*
*/		 13, 13,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14 \ /*
*/		 13, 13,  0,  1,  2,  3,  4,  0,  1,  2,  3,  4,  5, 6, 7, 8, 9)

mat SIMS = S[.,3...]
local simno = colsof(SIMS)


tempname ltenm1
postfile `ltenm1' xaxis time0_0 time0_1 time1_0 time1_1 /*
*/	x1 pr_s1_0_0 pr_s1_0_0_lo pr_s1_0_0_hi pr_s1_0_1 pr_s1_0_1_lo pr_s1_0_1_hi pr_s1_1_0 pr_s1_1_0_lo pr_s1_1_0_hi pr_s1_1_1 pr_s1_1_1_lo pr_s1_1_1_hi lte_s1 lte_s1_lo lte_s1_hi /*
*/	x2 pr_s2_0_0 pr_s2_0_0_lo pr_s2_0_0_hi pr_s2_0_1 pr_s2_0_1_lo pr_s2_0_1_hi pr_s2_1_0 pr_s2_1_0_lo pr_s2_1_0_hi pr_s2_1_1 pr_s2_1_1_lo pr_s2_1_1_hi lte_s2 lte_s2_lo lte_s2_hi /*
*/	x3 pr_s3_0_0 pr_s3_0_0_lo pr_s3_0_0_hi pr_s3_0_1 pr_s3_0_1_lo pr_s3_0_1_hi pr_s3_1_0 pr_s3_1_0_lo pr_s3_1_0_hi pr_s3_1_1 pr_s3_1_1_lo pr_s3_1_1_hi lte_s3 lte_s3_lo lte_s3_hi /*
*/	using "LTE Example--NM1.dta", replace


qui foreach c of numlist 1(1)`simno' {
	local x1 = S[1,`c']
	local x2 = S[2,`c']
	local x3 = S[3,`c']
	local t0_0 = S[4,`c']
	local t0_1 = S[5,`c']
	local t1_0 = S[6,`c']
	local t1_1 = S[7,`c']
	
	*** Y_{t-1} == 0 & Y_{t-1} == 1	
	foreach s of numlist 1(1)3 {
		foreach b of numlist 0 1 {
			foreach o of numlist 0 1 {
				tempvar prob_s`s'_`o'_`b'
				setx x `x`s'' t `t`o'_`b'' t2 `t`o'_`b''^2 t3 `t`o'_`b''^3
				simqi, prval(1) genpr(`prob_s`s'_`o'_`b'')
				sum `prob_s`s'_`o'_`b'', meanonly
				local pr_s`s'_`o'_`b' = r(mean)
				_pctile `prob_s`s'_`o'_`b'', perc(2.5 97.5)
				local pr_s`s'_`o'_`b'_lo = r(r1)
				local pr_s`s'_`o'_`b'_hi = r(r2)
			}
		}
	
		tempvar lteff_s`s'
		gen `lteff_s`s'' = `prob_s`s'_0_1' - `prob_s`s'_0_0'
		sum `lteff_s`s'', meanonly
		local lte_s`s' = r(mean)
		_pctile `lteff_s`s'', perc(2.5 97.5)
		local lte_s`s'_lo = r(r1)
		local lte_s`s'_hi = r(r2)
	}
	
	local xaxis = `c' - 2
	
	post `ltenm1' (`xaxis') (`t0_0') (`t0_1') (`t1_0') (`t1_1') /*
*/	(`x1') (`pr_s1_0_0') (`pr_s1_0_0_lo') (`pr_s1_0_0_hi') (`pr_s1_0_1') (`pr_s1_0_1_lo') (`pr_s1_0_1_hi') (`pr_s1_1_0') (`pr_s1_1_0_lo') (`pr_s1_1_0_hi') (`pr_s1_1_1') (`pr_s1_1_1_lo') (`pr_s1_1_1_hi') (`lte_s1') (`lte_s1_lo') (`lte_s1_hi') /*
*/	(`x2') (`pr_s2_0_0') (`pr_s2_0_0_lo') (`pr_s2_0_0_hi') (`pr_s2_0_1') (`pr_s2_0_1_lo') (`pr_s2_0_1_hi') (`pr_s2_1_0') (`pr_s2_1_0_lo') (`pr_s2_1_0_hi') (`pr_s2_1_1') (`pr_s2_1_1_lo') (`pr_s2_1_1_hi') (`lte_s2') (`lte_s2_lo') (`lte_s2_hi') /*
*/	(`x3') (`pr_s3_0_0') (`pr_s3_0_0_lo') (`pr_s3_0_0_hi') (`pr_s3_0_1') (`pr_s3_0_1_lo') (`pr_s3_0_1_hi') (`pr_s3_1_0') (`pr_s3_1_0_lo') (`pr_s3_1_0_hi') (`pr_s3_1_1') (`pr_s3_1_1_lo') (`pr_s3_1_1_hi') (`lte_s3') (`lte_s3_lo') (`lte_s3_hi')
}

postclose `ltenm1'

*****************************************************************
*** Variations of long-term effects based on the counterfactual of interest for non-monotonic 2 duration dependence
*****************************************************************
use "NM2.dta", clear

btscs y tim group, g(t) nspline(3)
replace t = (t + 1)
gen t2 = t^2
gen t3 = t^3

set seed 648
estsimp logit y x t t2 t3

************************ All Types of Effects *****************************
* Rows: x (1-time), x (temporary), x (permanent), t0_0, t0_1, t1_0, t1_1
* Columns: Baseline, counterfactual, then simulations
mat S = (0 , .5,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0 \ /*
*/  	 0 , .5, .5, .5, .5,  0,  0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0 \ /*
*/		 0 , .5, .5, .5, .5, .5, .5, .5, .5, .5, .5, .5,.5,.5,.5,.5,.5 \ /*
*/		 13, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28 \ /*
*/		 13, 13,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14 \ /*
*/		 13, 13,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14 \ /*
*/		 13, 13,  0,  1,  2,  3,  4,  0,  1,  2,  3,  4,  5, 6, 7, 8, 9)

mat SIMS = S[.,3...]
local simno = colsof(SIMS)


tempname ltenm2
postfile `ltenm2' xaxis time0_0 time0_1 time1_0 time1_1 /*
*/	x1 pr_s1_0_0 pr_s1_0_0_lo pr_s1_0_0_hi pr_s1_0_1 pr_s1_0_1_lo pr_s1_0_1_hi pr_s1_1_0 pr_s1_1_0_lo pr_s1_1_0_hi pr_s1_1_1 pr_s1_1_1_lo pr_s1_1_1_hi lte_s1 lte_s1_lo lte_s1_hi /*
*/	x2 pr_s2_0_0 pr_s2_0_0_lo pr_s2_0_0_hi pr_s2_0_1 pr_s2_0_1_lo pr_s2_0_1_hi pr_s2_1_0 pr_s2_1_0_lo pr_s2_1_0_hi pr_s2_1_1 pr_s2_1_1_lo pr_s2_1_1_hi lte_s2 lte_s2_lo lte_s2_hi /*
*/	x3 pr_s3_0_0 pr_s3_0_0_lo pr_s3_0_0_hi pr_s3_0_1 pr_s3_0_1_lo pr_s3_0_1_hi pr_s3_1_0 pr_s3_1_0_lo pr_s3_1_0_hi pr_s3_1_1 pr_s3_1_1_lo pr_s3_1_1_hi lte_s3 lte_s3_lo lte_s3_hi /*
*/	using "LTE Example--NM2.dta", replace


qui foreach c of numlist 1(1)`simno' {
	local x1 = S[1,`c']
	local x2 = S[2,`c']
	local x3 = S[3,`c']
	local t0_0 = S[4,`c']
	local t0_1 = S[5,`c']
	local t1_0 = S[6,`c']
	local t1_1 = S[7,`c']
	
	*** Y_{t-1} == 0 & Y_{t-1} == 1	
	foreach s of numlist 1(1)3 {
		foreach b of numlist 0 1 {
			foreach o of numlist 0 1 {
				tempvar prob_s`s'_`o'_`b'
				setx x `x`s'' t `t`o'_`b'' t2 `t`o'_`b''^2 t3 `t`o'_`b''^3
				simqi, prval(1) genpr(`prob_s`s'_`o'_`b'')
				sum `prob_s`s'_`o'_`b'', meanonly
				local pr_s`s'_`o'_`b' = r(mean)
				_pctile `prob_s`s'_`o'_`b'', perc(2.5 97.5)
				local pr_s`s'_`o'_`b'_lo = r(r1)
				local pr_s`s'_`o'_`b'_hi = r(r2)
			}
		}
	
		tempvar lteff_s`s'
		gen `lteff_s`s'' = `prob_s`s'_0_1' - `prob_s`s'_0_0'
		sum `lteff_s`s'', meanonly
		local lte_s`s' = r(mean)
		_pctile `lteff_s`s'', perc(2.5 97.5)
		local lte_s`s'_lo = r(r1)
		local lte_s`s'_hi = r(r2)
	}
	
	local xaxis = `c' - 2
	
	post `ltenm2' (`xaxis') (`t0_0') (`t0_1') (`t1_0') (`t1_1') /*
*/	(`x1') (`pr_s1_0_0') (`pr_s1_0_0_lo') (`pr_s1_0_0_hi') (`pr_s1_0_1') (`pr_s1_0_1_lo') (`pr_s1_0_1_hi') (`pr_s1_1_0') (`pr_s1_1_0_lo') (`pr_s1_1_0_hi') (`pr_s1_1_1') (`pr_s1_1_1_lo') (`pr_s1_1_1_hi') (`lte_s1') (`lte_s1_lo') (`lte_s1_hi') /*
*/	(`x2') (`pr_s2_0_0') (`pr_s2_0_0_lo') (`pr_s2_0_0_hi') (`pr_s2_0_1') (`pr_s2_0_1_lo') (`pr_s2_0_1_hi') (`pr_s2_1_0') (`pr_s2_1_0_lo') (`pr_s2_1_0_hi') (`pr_s2_1_1') (`pr_s2_1_1_lo') (`pr_s2_1_1_hi') (`lte_s2') (`lte_s2_lo') (`lte_s2_hi') /*
*/	(`x3') (`pr_s3_0_0') (`pr_s3_0_0_lo') (`pr_s3_0_0_hi') (`pr_s3_0_1') (`pr_s3_0_1_lo') (`pr_s3_0_1_hi') (`pr_s3_1_0') (`pr_s3_1_0_lo') (`pr_s3_1_0_hi') (`pr_s3_1_1') (`pr_s3_1_1_lo') (`pr_s3_1_1_hi') (`lte_s3') (`lte_s3_lo') (`lte_s3_hi')
}

postclose `ltenm2'

*********************************************************************
*** Average predictive differences
*********************************************************************

foreach s of numlist 1(1)4 {
	use "NPH Positive_APD_`s'.dta", clear

	gen scenario = `s'
	
	tempfile s`s'
	save `s`s'', replace
}

use `s1', clear

foreach s of numlist 2(1)4 {
	append using `s`s''
}

lab def scenario_l 1 "Beta = 0.2" 2 "Beta = 0.1" 3 "Beta = 0.04" 4 "Beta = 0.02"  
lab val scenario scenario_l

save "NwP--APD Summary.dta", replace


foreach s of numlist 1(1)4 {
	use "NPH Negative_APD_`s'.dta", clear

	gen scenario = `s'
	
	tempfile s`s'
	save `s`s'', replace
}

use `s1', clear

foreach s of numlist 2(1)4 {
	append using `s`s''
}


lab def scenario_l 1 "Beta = -0.2" 2 "Beta = -0.1" 3 "Beta = -0.04" 4 "Beta = -0.02"
lab val scenario scenario_l

save "PwN--APD Summary.dta", replace


****************************************************************
*** Predicted probability of dispute initiation for coalition governments and ideological distance
****************************************************************														
use "Clare.dta", clear

qui probit init_monad coalition dist_directional1 gov_ideology1 minority cap_1 prop_demborder1 prop_allyborder1 dep_trade veto_players peaceyears  peaceyears2 peaceyears3, robust
keep if e(sample)
sort ccode1 year

* How common are those observations that are statistically different from each other?
list ccode1 year quarter dist_directional1 if coalition == 1 & (inrange(dist_directional1, -125, -108) | inrange(dist_directional1, 92, 125))
tab ccode1 if coalition == 1 & (inrange(dist_directional1, -125, -75) | inrange(dist_directional1, 75, 125))

gen id = _n
sort id
tempfile c
save `c', replace

estsimp probit init_monad coalition dist_directional1 gov_ideology1 minority cap_1 prop_demborder1 prop_allyborder1 dep_trade veto_players peaceyears  peaceyears2 peaceyears3, robust

postfile pr_dist xaxis pr0 pr0_se pr0_lo pr0_hi pr1 pr1_se pr1_lo pr1_hi pr2 pr2_se pr2_lo pr2_hi using "Predicted Probabilities--Ideological Distance.dta", replace

qui foreach i of numlist -125(1)125 {
	setx mean
	setx gov_ideology1 0 dist_directional1 `i'
	
	simqi, prval(1) genpr(pr0) listx
	qui sum pr0
	local pr0 = r(mean)
	local pr0_se = r(sd)
	
	_pctile pr0, p(5.0 95.0)
	local pr0_lo = r(r1)
	local pr0_hi = r(r2)
	
	* Cohesive coalition
	setx mean
	setx gov_ideology1 0 dist_directional1 0
	
	simqi, prval(1) genpr(pr1) listx
	qui sum pr1
	local pr1 = r(mean)
	local pr1_se = r(sd)
	
	_pctile pr1, p(5.0 95.0)
	local pr1_lo = r(r1)
	local pr1_hi = r(r2)
	
	* Single-Party Majority Government
	setx mean
	setx gov_ideology1 0 dist_directional1 0 coalition 0 minority 0
	
	simqi, prval(1) genpr(pr2) listx
	qui sum pr2
	local pr2 = r(mean)
	local pr2_se = r(sd)
	
	_pctile pr2, p(5.0 95.0)
	local pr2_lo = r(r1)
	local pr2_hi = r(r2)
	
	drop pr0 pr1 pr2
	post pr_dist (`i') (`pr0') (`pr0_se') (`pr0_lo') (`pr0_hi') (`pr1') (`pr1_se') (`pr1_lo') (`pr1_hi') (`pr2') (`pr2_se') (`pr2_lo') (`pr2_hi')	
}

postclose pr_dist

use "Clare.dta", clear

qui probit init_monad coalition dist_directional1 gov_ideology1 minority cap_1 prop_demborder1 prop_allyborder1 dep_trade veto_players peaceyears  peaceyears2 peaceyears3, robust
keep if e(sample)
sort ccode1 year

gen id = _n
sort id

keep dist_directional1 peaceyears
save "dist.dta", replace


*************************************************************
*** Model fit for Clare (2010)
*************************************************************
use "Clare.dta", clear

qui probit init_monad coalition dist_directional1 gov_ideology1 minority cap_1 prop_demborder1 prop_allyborder1 dep_trade veto_players peaceyears  peaceyears2 peaceyears3, robust

*** Classification table
estat class

*** expected Percent Correctly Predicted (ePCP)
local dv `e(depvar)'

preserve
	keep if e(sample)
	local obs = _N
	predict pr, pr
	egen sumpr1 = sum(pr) if `dv' == 1
	egen sumpr0 = sum(1-pr) if `dv' == 0
	duplicates drop sumpr1 sumpr0, force
	egen ePCP_rtotal = rowtotal(sumpr1 sumpr0)
	egen ePCP_total = total(ePCP_rtotal)
	gen ePCP = 100 * (ePCP_total / `obs')
	duplicates drop ePCP, force
	qui sum ePCP
	local epcp = round(r(mean), 0.01)
	di "ePCP = `epcp'%"
restore	

*** Receiver-operating characteristics (ROC) curve and area under the curve (AUC)
qui probit init_monad coalition dist_directional1 gov_ideology1 minority cap_1 prop_demborder1 prop_allyborder1 dep_trade veto_players peaceyears  peaceyears2 peaceyears3, robust
predict p, pr

estat class

lroc
roctab init_monad p
local auc = round(r(area), 0.001)
local lo = round(r(lb), 0.001)
local hi = round(r(ub), 0.001)

cap postclose epcp
postfile epcp cut prop0 prop1 using "Expected Percent Correctly Predicted.dta", replace

qui sum p
local start = round(r(min), 0.001)
local end = round(r(max), 0.001)

qui foreach i of numlist `start'(0.001)`end' {
	estat class, cut(`i')
	local prop0 = 0.01 * r(P_n0)
	local prop1 = 0.01 * r(P_p1)
	
	post epcp (`i') (`prop0') (`prop1')
}

postclose epcp

preserve
	use "Expected Percent Correctly Predicted.dta", clear

	di "Area under the curve (AUC) = " `auc' "  95% CI = [" `lo' ", " `hi' "]"
	
	sort prop1
	twoway (line prop0 prop1) (function y=abs(x-1), lpattern(dash)), ytitle("Proportion of 0's Correctly Predicted") xtitle("Proportion of 1's Correctly Predicted") legend(off) 
restore

preserve
	use "Expected Percent Correctly Predicted.dta", clear
	
	local it = 1/(_N-1)
	egen x = fill(0(`it')1)
	gen y = abs(x-1)

	sort prop1
	save "ROC.dta", replace
restore

*** Separation plot (save the file here, do the rest in R)
preserve
	qui probit init_monad coalition dist_directional1 gov_ideology1 minority cap_1 prop_demborder1 prop_allyborder1 dep_trade veto_players peaceyears  peaceyears2 peaceyears3, robust
	keep if e(sample)
	keep init_monad p
	save "Separation Plot.dta", replace
restore


*****************************************************************
*** Variability of substantive effects
*****************************************************************
clear
set obs 1000

mat S = (0.0, 0.2, 0.4, 0.6, 0.8, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6 \ /*
*/		 0.2, 0.4, 0.6, 0.8, 1, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0)
local s = colsof(S)
local beta = 1

qui foreach i of numlist 1(1)`s' {
	local start = S[1,`i']
	local stop = S[2,`i']

	gen pr`i' = `start'+(`stop'-`start')*runiform()
	qui sum pr`i'
	local pr_mn = r(mean)
	local pr_sd = r(sd)
	local pr_min = r(min)
	local pr_max = r(max)
	
	gen me`i' = pr`i'*(1-pr`i')*`beta'
	qui sum me`i'
	local me_mn = r(mean)
	local me_sd = r(sd)
	local me_min = r(min)
	local me_max = r(max)	
}

gen id = _n
reshape long pr me, i(id) j(s)
gen s2 = s-5
sort s id

save "Variability.dta", replace




*************************************************************
*** Testing for Non-Proportional Hazards
*************************************************************
use "Clare.dta", clear

foreach v of varlist peaceyear* {
	gen dist_`v' = dist_directional1 * `v'
}

* Wald test
qui probit init_monad coalition dist_directional1 gov_ideology1 minority cap_1 prop_demborder1 prop_allyborder1 dep_trade veto_players peaceyears  peaceyears2 peaceyears3 dist_p*, robust
test dist_peaceyears = dist_peaceyears2 = dist_peaceyears3 = 0

* Likelihood-ratio test
qui probit init_monad coalition dist_directional1 gov_ideology1 minority cap_1 prop_demborder1 prop_allyborder1 dep_trade veto_players peaceyears  peaceyears2 peaceyears3 dist_p*
estimates store full_model

qui probit init_monad coalition dist_directional1 gov_ideology1 minority cap_1 prop_demborder1 prop_allyborder1 dep_trade veto_players peaceyears  peaceyears2 peaceyears3
estimates store c_model
lrtest full_model c_model

